In [1]:
# Import necessary libraries
import pandas as pd
import numpy as np
import math
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.interpolate import interp1d
import yaml
import os
import scanpy as sc
import plotly.express as px

# Load the imaging data

Load imaging data¶

In [2]:
figDir = "./figures"
if not os.path.exists(figDir):
    os.makedirs(figDir)
outdir = "../data/output"
# Load colors dict
with open("../data/resources/iPSC_lines_map.yaml", 'r') as f:
    iPSC_lines_map = yaml.load(f, Loader=yaml.FullLoader)["lines"]
colorsmap = dict(zip([i["newName"] for i in iPSC_lines_map.values()],[i["color"] for i in iPSC_lines_map.values()]))
In [3]:
dataImaging=pd.read_csv('../data/imaging_data/organoidMultiplexing_growthCurves_quant.csv')

# Filter and reformat the imaging data
dataImaging = dataImaging.loc[dataImaging["FileName"].str.contains("MIX"), ["Line","Day","LineRep","EquivalentDiameterArea(microns)"]]
dataImaging["replicate"] = dataImaging["LineRep"].str.split("_",expand=True)[1].astype(int)
dataImaging = dataImaging[["Line","Day","replicate","EquivalentDiameterArea(microns)"]].rename(columns={"Line":"Mix"})
# removing incomplete mix 7
dataImaging = dataImaging[dataImaging["Mix"] != "MIX7"]

Load census data¶

In [4]:
dataCensus=pd.read_csv('../data/census_seq/CensusSeq_data.csv')
# Filter and reformat the census data
dataCensus["Timepoint"] = dataCensus["Timepoint"].str.split(" ", expand= True)[1].astype(int)
dataCensus = dataCensus[["Sample name","DONOR","REPRESENTATION","MIX ID","Timepoint"]].rename(columns={"MIX ID":"Mix","Timepoint":"Day"})
dataCensus["Mix"] = "MIX"+dataCensus["Mix"].astype(str)
# removing incomplete mix 7
dataCensus = dataCensus[dataCensus["Mix"] != "MIX7"]
dataCensus["DONOR"] = dataCensus["DONOR"].replace({"CTL01A":"CTL01"})
dataCensus["DONOR"] = dataCensus["DONOR"].replace({"CHD2WT":"UCSFi001-A"})
dataCensus["DONOR"] = dataCensus["DONOR"].replace({"CHD8WT":"H9"})
In [5]:
# Define the density volume
densityVolume = 0.000767495

# Transform diameter to cells
dataImaging["EquivalentSphereVolume"] = (np.power((dataImaging["EquivalentDiameterArea(microns)"]/2), 3))*math.pi*(4/3)
dataImaging["nCells"] = dataImaging["EquivalentSphereVolume"]*densityVolume
#%%
# Compute DeltaGrowth
GrowthRate = pd.DataFrame()
# Loop over each unique mix in the census data
for m in dataCensus.Mix.unique().tolist():
    GrowthRateMix = pd.DataFrame()
    # Filter the data for the current mix
    DataMix = dataCensus[dataCensus["Mix"] == m]
    DataMix = DataMix[DataMix["Day"].isin([5,12,-2])]
    DataMix["ClosestDay"] = DataMix["Day"].replace({5:4,12:10,-2:0})
    TPTS = np.sort(DataMix.ClosestDay.unique())
    # Loop over each unique timepoint for the current mix
    for t in TPTS:
        # Compute the average representation for each donor at the current timepoint
        DataMixTPT = DataMix[DataMix["ClosestDay"] == t].groupby("DONOR")["REPRESENTATION"].mean().to_frame().T.reset_index(drop=True)
        DataMixTPT.columns = ["AvgCensus."+i for i in DataMixTPT.columns.tolist()]
        # Compute the mean and standard deviation of the number of cells for the current mix and timepoint
        Mean = round(dataImaging.loc[(dataImaging["Day"] == t) &(dataImaging["Mix"] == m),"nCells"].mean(),2)
        std = round(dataImaging.loc[(dataImaging["Day"] == t) &(dataImaging["Mix"] == m),"nCells"].std(),2)
        # Create a dataframe with the computed values
        NcellsDF = pd.DataFrame({"Mix":m,"Day":t,"MeanCells":Mean, "STDcells":std}, index=[0])
        NcellsDF = pd.concat([NcellsDF,DataMixTPT], axis = 1)
        # Append the dataframe to the growth rate dataframe for the current mix
        GrowthRateMix = pd.concat([GrowthRateMix,NcellsDF], ignore_index=True)
    # If there are at least two timepoints for the current mix, compute the growth rate
    if GrowthRateMix.shape[0] >= 2:
        # Compute the growth rate and add it to the growth rate dataframe
        GrowthRateMix = GrowthRateMix.sort_values("Day")
        GrowThaRateVals = GrowthRateMix.T.loc[DataMixTPT.columns]
        GrowThaRateVals = GrowThaRateVals.T.div(GrowThaRateVals.T.shift(1)).T.dropna(axis = 1)
        GrowThaRateVals.index = ["GR."+i.split(".")[1] for i in GrowThaRateVals.index]
        GrowThaRateVals.columns = ["D{}_to_D{}_growthRate".format(GrowthRateMix.loc[i,"Day"],GrowthRateMix.loc[i+1,"Day"]) for i in GrowthRateMix.sort_values("Day").index.tolist()[:-1]]
        GrowThaRateVals["Mix"] = m
        GrowThaRateVals["DONOR"] = [i.split(".")[1] for i in GrowThaRateVals.index.tolist()]
        # Add the mean and standard deviation of the number of cells and the average census to the growth rate dataframe
        for d in GrowthRateMix.Day.unique():
            GrowThaRateVals["MeanCells.D{}".format(d)] = GrowthRateMix.loc[GrowthRateMix["Day"] == d,"MeanCells"].values[0]
            GrowThaRateVals["STDcells.D{}".format(d)] = GrowthRateMix.loc[GrowthRateMix["Day"] == d,"STDcells"].values[0]
            AvgCensus = GrowthRateMix.loc[GrowthRateMix["Day"] == d, [i for i in GrowthRateMix.columns.tolist() if i.startswith("Avg")]].T
            AvgCensus.index = [i.split(".")[1] for i in AvgCensus.index.tolist()]
            GrowThaRateVals["AvgCensus.D{}".format(d)] = AvgCensus.loc[GrowThaRateVals.DONOR].to_numpy().flatten()
        # Append the growth rate dataframe for the current mix to the overall growth rate dataframe
        GrowthRate = pd.concat([GrowthRate,GrowThaRateVals], ignore_index=True)

D0 to D4 growth rate stability¶

In [6]:
fig, axes = plt.subplots(1,1, figsize=(15,5), dpi=300, sharex=False,sharey=False)
sns.boxplot(data=GrowthRate, x="DONOR", y="D0_to_D4_growthRate", whis=(0, 100),palette=colorsmap)
sns.stripplot(data=GrowthRate, x="DONOR", y="D0_to_D4_growthRate", size=10, color=".3")
axes.title.set_text("Growth rate stability from D0 to D4")
axes.set_xlabel("Cell line", fontsize=20)
axes.set_ylabel("D0 to D4 growth rate", fontsize=20)
sns.despine(offset=30)
plt.show()

fig.savefig(figDir+"/GrowthRateStability.0to4.svg", format='svg', bbox_inches='tight')

D4 to D10 growth rate stability¶

In [7]:
fig, axes = plt.subplots(1,1, figsize=(15,5), dpi=300, sharex=False,sharey=False)

sns.boxplot(data=GrowthRate, x="DONOR", y="D4_to_D10_growthRate", whis=(0, 100),palette=colorsmap)
sns.stripplot(data=GrowthRate, x="DONOR", y="D4_to_D10_growthRate", size=10, color=".3")
axes.title.set_text("Growth rate stability from D4 to D10")
axes.set_xlabel("Cell line", fontsize=20)
axes.set_ylabel("D4 to D10 growth rate", fontsize=20)
sns.despine(offset=30)
plt.show()
fig.savefig(figDir+"/GrowthRateStability.4to10.svg", format='svg', bbox_inches='tight')
In [8]:
# %%

sns.histplot(GrowthRate["D4_to_D10_growthRate"])
total=GrowthRate["D4_to_D10_growthRate"].tolist()
total.extend(GrowthRate["D0_to_D4_growthRate"].tolist())
sns.histplot(total)
Out[8]:
<AxesSubplot: xlabel='D4_to_D10_growthRate', ylabel='Count'>
In [9]:
import matplotlib.pyplot as plt

# Assuming GrowthRate is a pandas DataFrame
D4_to_D10_growthRate = GrowthRate["D4_to_D10_growthRate"]
D0_to_D4_growthRate = GrowthRate["D0_to_D4_growthRate"]

plt.hist([D4_to_D10_growthRate, D0_to_D4_growthRate,total], label=['D4_to_D10_growthRate', 'D0_to_D4_growthRate', 'Total'], stacked=False,density=True)
plt.legend(loc='upper right')
plt.show()
In [10]:
sns.histplot(total, stat='probability')
# Assume 'ax' is the axes object containing the histogram
ax = plt.gca()

# Get the rectangles representing the bars
patches = ax.patches

# Extract the heights (frequencies) and the left edges of the rectangles
frequencies = [patch.get_height() for patch in patches]
bin_edges = [patch.get_xy()[0] for patch in patches]

# Calculate the probabilities
probabilities = frequencies
In [11]:
# %%
mix_number=10
cell_dict={}
EB_cell_number=20000
total_cycles=10000
for i in range(1,total_cycles):
    generated_data = np.random.choice(bin_edges, size=mix_number, p=probabilities)
    cell_dict[i]=EB_cell_number*generated_data

simulated_exp=pd.DataFrame(cell_dict)

simulated_exp = simulated_exp.div(simulated_exp.sum(axis=0), axis=1)>0.05
sns.histplot(simulated_exp.sum(axis=0))
Out[11]:
<AxesSubplot: ylabel='Count'>
In [12]:
min_x = 0
max_x = 6

sns.histplot(total, stat='probability')
# Assume 'ax' is the axes object containing the histogram
ax = plt.gca()

# Get the rectangles representing the bars
patches = ax.patches

# Extract the heights (frequencies) and the left edges of the rectangles
frequencies = [patch.get_height() for patch in patches]
bin_edges = [patch.get_xy()[0] for patch in patches]

# Calculate the probabilities
probabilities = frequencies

# Get density from Seaborn
x, y = sns.kdeplot(np.random.choice(bin_edges, size=10000, p=probabilities),clip=(min_x, max_x)).lines[0].get_data()

probabilities=y/sum(y)

mix_number=10
cell_dict={}
EB_cell_number=20000
total_cycles=10000
for i in range(1,total_cycles):
    generated_data = np.random.choice(x, size=mix_number, p=probabilities)
    cell_dict[i]=EB_cell_number/mix_number*generated_data

simulated_exp=pd.DataFrame(cell_dict)

simulated_exp = simulated_exp.div(simulated_exp.sum(axis=0), axis=1)>0.05
sns.histplot(simulated_exp.sum(axis=0),binwidth=1)
Out[12]:
<AxesSubplot: ylabel='Probability'>
In [13]:
def empirical_dist(total, min_x, max_x):
    plt.clf()
    sns.histplot(total, stat='probability')
    ax = plt.gca()
    patches = ax.patches
    frequencies = [patch.get_height() for patch in patches]
    bin_edges = [patch.get_xy()[0] for patch in patches]
    x, y = sns.kdeplot(np.random.choice(bin_edges, size=10000, p=frequencies),clip=(min_x, max_x)).lines[0].get_data()
    probabilities=y/sum(y)
    return x, probabilities

cell_dict={}
median_dict={}


min_x = 0
max_x = 9

mix_number=2
max_mix=30
EB_cell_number=20000
total_cycles=10000
nFinalCells = 100000
minCellsPerCelltype = 100
RarestCelltypeRate = 0.10
#what is the minimum number of cells per genotype given RarestCelltypeRate and  minCellsPerCelltype required?
# RarestCelltypeRate(x100) : 100 = minCellsPerCelltype : X
MinCellsPerGenotype = (100*minCellsPerCelltype)/(RarestCelltypeRate*100)
MinCellsPerGenotype_Rate = MinCellsPerGenotype/nFinalCells



x, probabilities = empirical_dist(GrowthRate["D0_to_D4_growthRate"].tolist(), min_x, max_x)
x1, probabilities1 = empirical_dist(GrowthRate["D4_to_D10_growthRate"].tolist(), min_x, max_x)

k=0
while mix_number <= max_mix:
    cell_dict[mix_number]={}
    for i in range(1,total_cycles):
        generated_data = np.random.choice(x, size=mix_number, p=probabilities)
        generated_data_2 = np.random.choice(x1, size=mix_number, p=probabilities1)
        if sum(generated_data)>mix_number:
            k+=1
            generated_cell=(EB_cell_number/mix_number)*generated_data
            # if sum(generated_cell*generated_data_2)>sum(generated_cell):
            #     cell_dict[mix_number][k]=generated_cell*generated_data_2
            cell_dict[mix_number][k]=generated_cell
    

    simulated_exp=pd.DataFrame(cell_dict[mix_number])


    plt.clf()
    plot_data=(simulated_exp.div(simulated_exp.sum(axis=0), axis=1)>MinCellsPerGenotype_Rate).sum(axis=0)
    testplot=sns.histplot(plot_data,discrete=True,stat='probability')
    testplot.set(xlabel='Number of analysable genotypes', ylabel='Frequency')
    testplot.set(xticks=range(0,15,1))
    plt.title('Number of genotypes mixed: '+str(mix_number))
    # plt.show()
    plt.savefig(figDir+'/histogram_dual_'+str(mix_number)+'.png')
    
    median_dict[mix_number]=[plot_data.mean()-plot_data.std(),plot_data.mean(),plot_data.mean()+plot_data.std()]
    mix_number+=1
    print(plot_data.std())

plt.clf()



data_median=pd.DataFrame(median_dict).T
data_median['x']=data_median.index
0.2720136693772788
0.4014987651198087
0.5142667238637149
0.6290950269315375
0.7497535791375715
0.8454721343492219
0.9471576840718762
1.0267961239790546
1.1008441298980556
1.1862231578086202
1.262627302682708
1.3242728947384765
1.3884062866579954
1.4222098214165106
1.4744492139393681
1.5347643812676524
1.5789831168366635
1.6260193020108062
1.6899430581407442
1.697537004618319
1.7624928630124848
1.809411680664435
1.8540251868661528
1.841915994137707
1.9218028072236684
1.9281266789300264
1.9530723989390903
2.017881002712317
2.0282749255674513
<Figure size 432x288 with 0 Axes>
In [14]:
#%%
# PLot empirical distribution
sns.set_style("whitegrid")
sns.set_style("ticks")
fig,ax = plt.subplots(ncols=1,nrows=1, figsize=(7, 5),dpi=100 )

sns.set_style("white")
#plt.plot(data_median["x"], data_median[1], '-', linewidth = 3, label = "Median")
sns.lineplot(data=data_median,x="x",y=1,color="#FDA63A",
                 linewidth=3, ax = ax)
sns.lineplot(data=data_median,x="x",y=0,color="blue",
                 linewidth=1.5, ax = ax)
sns.lineplot(data=data_median,x="x",y=2,color="blue",
                 linewidth=1.5, ax = ax)
plt.fill_between(data_median["x"], data_median[0], data_median[2] , 
                 alpha=0.2, color='blue',linewidth=0.0)
#plt.xlabel('pseudotime', size=30)

ax.set_ylim([0,
             round(data_median[2].max()+1) ])

ax.set_xlim([0, data_median["x"].max()])

ax.set_xticks([0, 2,5,10,20])

ax.set_xlabel('Mixed cell lines') # Y label
ax.set_ylabel('Recovered cell lines') # Y label

sns.despine(fig=ax.figure)
fig.savefig(figDir+"/MinGenotypeRate0.05.Simulation.svg", format='svg', bbox_inches='tight')

Assess main cell types abundances¶

In [15]:
adataComp = sc.read(outdir+"/adatas/adataPaga.h5ad")
In [16]:
leidenOrder=["ProliferatingProgenitors","RadialGliaProgenitors","OuterRadialGliaAstrocytes","CajalR_like","Neurons","MigratingNeurons","GlutamatergicNeurons_early","GlutamatergicNeurons_late","Interneurons_GAD2","Interneurons"]

for paradigm in adataComp.obs.type.unique():
    adataCompStage = adataComp[adataComp.obs["type"] == paradigm]
    compositions = pd.DataFrame(adataCompStage.obs.groupby(["leidenAnnotated","stage"]).size())
    compositions = compositions.reset_index().rename(columns={0:"number_of_cells"})
    compositions["cells_fraction"] = compositions["number_of_cells"] / np.array(compositions.groupby("stage")["number_of_cells"].sum()[compositions["stage"].tolist()])
    compositions["cells_fraction"] = round(compositions["cells_fraction"],2)
    fig = px.bar(compositions, x="stage", y="cells_fraction", color="leidenAnnotated", hover_data=compositions, text="cells_fraction",
                 category_orders={"stage":['early', 'mid', 'late'],
                                 "leidenAnnotated":leidenOrder}, height=800,width=1000, template="plotly_white",title=paradigm,
                 color_discrete_map = dict(zip(adataCompStage.obs["leidenAnnotated"].cat.categories,adataCompStage.uns["leidenAnnotated_colors"]))
    )

    fig.update_traces(textangle=0, marker_line_color='rgb(8,48,107)', textposition = 'inside',
                    marker_line_width=1.5, opacity=1)



    fig.update_layout(font=dict(size=25, color='black'),
        yaxis = dict(tickfont = dict(size=30)),
        xaxis = dict(tickfont = dict(size=30)))

    fig.show()

~ 10 % rarest cell type should be enough to capture main cell types¶

In [17]:
cell_dict={}
median_dict={}


min_x = 0
max_x = 9

mix_number=2
max_mix=20
EB_cell_number=20000
total_cycles=10000
nFinalCells = 100000
minCellsPerCelltype = 100
RarestCelltypeRate = 0.09
#what is the minimum number of cells per genotype given RarestCelltypeRate and  minCellsPerCelltype required?
# RarestCelltypeRate(x100) : 100 = minCellsPerCelltype : X
MinCellsPerGenotype = (100*minCellsPerCelltype)/(RarestCelltypeRate*100)
MinCellsPerGenotype_Rate = MinCellsPerGenotype/nFinalCells
print(MinCellsPerGenotype_Rate)


x, probabilities = empirical_dist(GrowthRate["D0_to_D4_growthRate"].tolist(), min_x, max_x)
x1, probabilities1 = empirical_dist(GrowthRate["D4_to_D10_growthRate"].tolist(), min_x, max_x)

k=0
while mix_number <= max_mix:
    cell_dict[mix_number]={}
    for i in range(1,total_cycles):
        generated_data = np.random.choice(x, size=mix_number, p=probabilities)
        generated_data_2 = np.random.choice(x1, size=mix_number, p=probabilities1)
        if sum(generated_data)>mix_number:
            k+=1
            generated_cell=(EB_cell_number/mix_number)*generated_data
            # if sum(generated_cell*generated_data_2)>sum(generated_cell):
            #     cell_dict[mix_number][k]=generated_cell*generated_data_2
            cell_dict[mix_number][k]=generated_cell
    

    simulated_exp=pd.DataFrame(cell_dict[mix_number])


    plt.clf()
    plot_data=(simulated_exp.div(simulated_exp.sum(axis=0), axis=1)>MinCellsPerGenotype_Rate).sum(axis=0)
    testplot=sns.histplot(plot_data,discrete=True,stat='probability')
    testplot.set(xlabel='Number of analysable genotypes', ylabel='Frequency')
    testplot.set(xticks=range(0,15,1))
    plt.title('Number of genotypes mixed: '+str(mix_number))
    # plt.show()
    
    median_dict[mix_number]=[plot_data.mean()-plot_data.std(),plot_data.mean(),plot_data.mean()+plot_data.std()]
    mix_number+=1
    print(plot_data.std())

plt.clf()



data_median=pd.DataFrame(median_dict).T
data_median['x']=data_median.index


#%%
# PLot empirical distribution
sns.set_style("whitegrid")
sns.set_style("ticks")
fig,ax = plt.subplots(ncols=1,nrows=1, figsize=(7, 5),dpi=100 )

sns.set_style("white")
#plt.plot(data_median["x"], data_median[1], '-', linewidth = 3, label = "Median")
sns.lineplot(data=data_median,x="x",y=1,color="#FDA63A",
                 linewidth=3, ax = ax)
sns.lineplot(data=data_median,x="x",y=0,color="blue",
                 linewidth=1.5, ax = ax)
sns.lineplot(data=data_median,x="x",y=2,color="blue",
                 linewidth=1.5, ax = ax)
plt.fill_between(data_median["x"], data_median[0], data_median[2] , 
                 alpha=0.2, color='blue',linewidth=0.0)
#plt.xlabel('pseudotime', size=30)

ax.set_ylim([0,
             round(data_median[2].max()+1) ])

ax.set_xlim([0, data_median["x"].max()])

ax.set_xticks([0, 2,5,10,20])

ax.set_xlabel('Mixed cell lines') # Y label
ax.set_ylabel('Recovered cell lines') # Y label

sns.despine(fig=ax.figure)
fig.savefig(figDir+"/MinGenotypeRate0.05.Simulation.RarestCT:{}.svg".format(RarestCelltypeRate), format='svg', bbox_inches='tight')
0.011111111111111112
0.2793244122848294
0.4127758222493213
0.54609069067006
0.6512412831191534
0.7809656090841938
0.8702182134538637
0.9618361064033875
1.0509956284952913
1.137219942219343
1.2135456161374711
1.2609719963584933
1.3181411779107048
1.3764093894239915
1.4590027628053408
1.4844533129696897
1.54647433569609
1.6124101501398906
1.6141509942954693
1.6851675743182277
<Figure size 432x288 with 0 Axes>
In [18]:
from matplotlib import pylab

from matplotlib.pyplot import figure
sns.set_style("whitegrid")
pylab.rcParams['figure.figsize'] = (5, 7)

FittedValues["MixedGenotypes"] = FittedValues.index.tolist()
plot = FittedValues.melt("MixedGenotypes", value_name="Recovered Genotypes").copy()
plot = plot[plot["variable"] == "Fitted_log_RarestCelltypeRate:{}".format(RarestCelltypeRate)]

sns.set_style("ticks")
fig,ax = plt.subplots(ncols=1,nrows=1, figsize=(7, 5),dpi=100 )
sns.lineplot(data=plot,x="MixedGenotypes",y="Recovered Genotypes",hue="variable",palette="Set2",
                 linewidth=3, ax = ax)

ax.get_legend().remove()# To remove the legend

# ax.set_ylabel('SHH expressing cells', fontsize = 30) # Y label
# ax.set_xlabel('genotype', fontsize = 30) # X label
ax.set_ylim([0,
             round(plot["Recovered Genotypes"].max()+1) ])

ax.set_xlim([0, (max_mix+1)])


#ax.set_yticks([0,.25,.5,.75,1])
fig.suptitle('')
ax.set_xlabel('Mixed cell lines') # Y label
ax.set_ylabel('Recovered cell lines') # Y label
ax.set_xticks([0, 2,5,10,20])



### Pick the vertical lines to plot
vlinelist = [5,10,20]
#fittedParams = curve_fit(powerFun, finalDF.index.tolist(),finalDF["Fitted_log_RarestCelltypeRate:{}".format(RarestCelltypeRate)])[0]
predRecover = powerFun(np.array(vlinelist), *fittedParams)
vlineDict = dict(zip(vlinelist,predRecover))
vlineDictDF = pd.DataFrame(vlineDict, index=["Recovered"]).T
vlineDictDF["mixedLines"] = vlineDictDF.index

for vl in list(vlineDict.keys()):
    plt.vlines(x=vl,ymax=vlineDict[vl],ymin=0, color='#e6e6e6',linestyle='--')
    plt.hlines(y=vlineDict[vl],xmax=vl,xmin=0, color='magenta',linestyle='-',linewidth = .1 )

ax.set_title('Number of recovered cell lines with more than {}% representation'.format(MinCellsPerGenotype_Rate*100))


sns.scatterplot(data = vlineDictDF, x = "mixedLines", y = "Recovered", color = "magenta", 
                s = 50,ax = ax, zorder=100)
ax.get_legend().remove() 

sns.despine(fig=ax.figure)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-18-77f59d0bc0dc> in <module>
      5 pylab.rcParams['figure.figsize'] = (5, 7)
      6 
----> 7 FittedValues["MixedGenotypes"] = FittedValues.index.tolist()
      8 plot = FittedValues.melt("MixedGenotypes", value_name="Recovered Genotypes").copy()
      9 plot = plot[plot["variable"] == "Fitted_log_RarestCelltypeRate:{}".format(RarestCelltypeRate)]

NameError: name 'FittedValues' is not defined